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ABSTRACT 

The  bounds  of  Boyd  J[1959)  and  Soms  ( 1980 1984)  for  the  tail  areas  of  the 
normal  and  t-distributions  are  used  to  obtain  a  new  method  of  approximating 
the  tail  areas.  The  absolute  and  relative  errors  are  given. 
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SIGNIFICANCE  AND  EXPLANATION 

A  reasonably  accurate  and  simple  method  of  evaluating  normal  and  t-tail 
areas  and  percentiles  is  described.  This  approach  controls  both  absolute  and 
relative  errors,  unlike  the  usual  methods  which  only  control  the  absolute 
error.  The  control  of  both  errors  is  of  importance  in  calculating  Bonferroni 
percentiles  and  Bonferroni  descriptive  probability  levels,  which  are  useful  in 
simultaneous  inference.  For  example,  a  nonspecialist  in  statistics,  unaware 
of  the  limitations  of  standard  routines,  might  report  a  Bonferroni  descriptive 
probability  level  that  had  a  relative  error  of  100%  or  more. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MR C,  and  not  with  the  author  of  this  report. 
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A  NEW  METHOD  OK  APPROXIMATING  NORMAL  AND  t-TAIL  AREAS 

Andrew  P.  Soms* 

1 .  THE  METHOD 

We  begin  by  introducing  notation  and  stating  the  main  results 
of  Boyd  (1959)  and  Soms  (1980,  1984).  Boyd  (1959)  showed  that  if 


_  ^  2  — 

<M  x )  ( 2n  )  exp(-x  /2),  F(x)  “  t  A(t)dt  , 


and  Rx  =  F(x)/'f’(x),  x  >  0,  then 


p (  x ,  Y  )  <  R  <  p(  x ,  Y  )  , 

*  '  min  x  K  max  ' 

where  p(x,Y)  =  (Y  +  1)/[(x2  +  ( 2/n ) (Y  +  I)2)  1/2  +  Yx] , 

Y  =  2/(ir  -  2),  Y.  =  it  -  1 ,  and  the  bounds  are  the  best 
max  min 

possible  in  the  class  fp(x,Y),  Y  >  -iK  This  is  also  discussed 
in  Johnson  and  Kotz  (1970,  Ch.  33).  Soms  (1980),  1984)  extended 
the  above  results  and  showed  that  if  for  arbitrary  real  k  >  0 
and  x  >  0, 

f  (t)  «  CM  ♦  t2/k)"(k+1>/2,  ck  =  T1SJL  +  UJJLL  , 

r(k/2)(irk)  ' 


F,  (x)  =  1  -  F  (X)  =  f  f  (t)dt  , 
k  k  k 

x 


iyx)  =  F  <  x )  /  [  (  1  +  xZA)f  (X)]  , 


2  2  k 

for  k  >  2,  Y  =  4c,  /(  1  -  4c.)  and  Y  .  =  - 

'  max  k  k  min 


-  1, 


2 ( k  +  2 )c. 


and  for  k  <  2,  Ym^n  and  Ymax  are  interchan<?ed»  and 
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p(  x  ,  Y  )  = 


_ 1  -t  V _ 

,2  2,  ,  2  1/2 

( x  Me  (1  t  y)  )  +  yx 

K 


Mi 


P ( X , Y  .  )  <  R  (x)  <  p( x , Y  )  , 

nun  k  max 


or  equivalently, 


-  •  , 


f  1  +  ~^f.  (x)P(x,Y  )  <  F  (x)  <  f  1  +  -p— If,  (x)p(x,Y  )  , 
k  k  min  k  k  k  max 

and  the  hounds  aqain  are  best  in  the  same  sense  as  for  the  normal. 


It  was  also  shown  there  that  if  k  =  2,  y  =  y  . 

max  mm 


Y  2  and 


Rk(x)  =  p(x,Y„). 

The  numerical  properties  of  these  bounds  are  discussed  in  the 
above  references.  The  important  fact  to  be  noted  here  is  that  the 
bounds  control  both  absolute  and  relative  error.  Usinq  the  hounds 
as  a  start inq  point  we  now  develop  a  simple  method  of  evaluatinq 
normal  and  t-tail  areas  that  controls  both  absolute  and  relative 
error,  as  opposed  to  the  usual  methods,  which  generally  only 
control  absolute  error. 

We  consider  estimates  of  the  tail  area  of  the  form 

fa— ;  jXV(*,y  ■  )Mx)  +  fl  -  a  k*fp(x,Y  )A(x)  (1.1) 

c  +  dx  min  c  +  dx  c  max 


for  the  tail  area  of  the  normal  and 


r  a  +  bx 


-^p(  x,  Y  .  )g  (x)  +  fl  -  3  *xlp(x,Y  )g  (x)  (1.2) 

+  dx  min  k  c  +  dx  max  k 


for  the  tail  area  of  the  t.,  where  cg-^.  ( >c )  -  (1  +  x^/k)fk(x).  vie 
want  the  estimates  to  lie  between  the  upper  and  lower  bounds  for 
the  tail  area  and  be  strictly  decreasing  functions  of  x  and 
therefore  impose  the  added  restrictions  that 

be  >  ad 


.  •  .*• 


o  <  c  1  , 

c  +  dx 

all  x  >  0.  Since  f ( 0 )  =  a/c,  we  may,  without  loss  of 
generality,  assume  that  c  =  1  and  so  our  weight  functions  f 
are  of  the  type 


f(x)  = 


a _ + _ bx 

n  dx  ' 


(1.3) 


where  0  <  a  <  1,  d>0,  b  >  ad  and  h/d  <1.  We  then  seek 

t-hat  particular  choice  of  f  which  minimizes  the  absolute 
error.  A  direct  computer  search  led  to 


.7  lx 

1  +  .7  lx 


(1.4) 


for  the  normal  and 


f  ( x )  = 


1  +  b  x  ' 
k 


(1.5) 


.71  +  l.fl/k  -  . 19/k^  , 


(1.6) 


for  Hie  t,  where,  as  noted  before,  k  is  the  degrees  of 
freedom.  (1.6)  was  obtained  by  finding  the  optimal  constants 
for  k  =  25,  10,  5,  3,  1.5,  1,  .5  and  fitting  a  constrained 
(constant  equal  to  .71)  regression  line.  However,  in  the 
interests  of  simplicity,  for  k  <  2,  we  did  not  interchange 

y  and  y  and  so  (1.5)  and  (1.6)  are  understood  to  apply 

mi n  max 


for  all  k  with  y  .  and  y 

min  max 


defined  as  for  k  >  2. 


Numerical  evidence  indicates  that,  at  least  for  k  =  1,  the  above 
optimal  estimate  is  still  a  decreasing  function  of  x. 

The  maximum  absolute  and  relative  errors  of  the  optimal 
estimates  are  remarkably  constant  over  the  range  1  <  k  <  00  and 
hence  we  only  oive  the  normal  figures.  For  (1.4),  the  maximum 
absolute  error  is  .66  x  10”^  and  the  maximum  relative  error  is 
.97  x  10_3.  we  emphasize  once  more,  that,  unlike  the  usual 
methods,  which  generally  control  only  absolute  error  and  give 
misleading  results  for  large  x  or  small  tail  areas,  the  above 


controls  both  absolute  and  relative  error  and  hence  can  be  used  to 
calculate  ordinary  and  Bonferroni  descriptive  levels  and  ordinary 
and  Bonferroni  percentiles. 

As  a  check,  we  calculated  the  standard  textbook  table  of  the 
normal,  given,  e.g.,  in  Brown  and  Hollander  (1977)  and  found  at 
most  a  difference  of  1  in  the  fourth  decimal  place.  We  also 
compared  the  small  normal  percentiles  given  in  Abramowitz  and 
Stegun  (1q65,  p.  977)  to  the  ones  obtained  from  (1.4)  and  after 
roundinu  both  to  three  decimal  places  found  that  there  was  at  most 
a  difference  of  1  in  the  third  decimal  place.  Similar  results 
apply  to  the  t. 

2.  CONCLUDING  RFMARXS 

We  have  given  a  method  of  calculating  normal  and  t-tail  areas 
and  percentiles  which  controls  both  absolute  and  relative  errors 
and  is  simple  to  implement  even  on  small  computers.  Four  short 
FORTRAN  programs  that  do  this  are  given  in  the  Appendix.  The 
programs  are  completely  self-contained,  using  a  short  FORTRAN 
subroutine  to  calculate  ck  based  on  Stirling’s  formula  for 
large  k  and  an  approximation  to  the  gamma  function  given  in 
Abramowitz  and  Stegun  (1965,  p.  257)  for  small  k. 


APPENDIX 


The  list.inqs  of  the  four  programs  referred  to  in  the  paper 
are  given  below.  If  k  <  2  and  n  is  small,  the  t-percentile 
program  may  be  slow. 

SUBROUTINE  NORT(X,TP) 

C  NORMAI,  TAI I.  AREA  FROM  X  TO  INFINITY,  X  NON- N  EC  AT  I V  E 
DATA  PIE/3. 1415427/ 

PHI( X)  =  1 . /SQRT( 2  *  P I E ) *  EX  P ( -X**2/2 . ) 

P(X,G)  =(  1  .  +  C.)/(  (X**2+2./PIF,*(  1  .+G)**2)**.5+G*X) 

PAMIN (X)=.71*X/( 1.+.71*X) 

TPL(X,GL,GU)=(PAMIN{X)*P(X,CL)+( 1 .-PAMIN ( X ) ) *P( X ,GU ) )*PHI(X) 
GMIN=PIE-1 . 

GMAX=2. /(PIE-2. ) 

TP=TPL( X ,GMIM ,GMAX ) 

RETURN 

END 

SUBROUTINE  NORTI ( ALPHA, XALP) 

C  UPPER  NORMAL  PERCENTILE  XALP  CORRESPONDING  TO  TAIL 
C  PROBABILITY  ALPHA,  ALPHA  BETWEEN  0  AND  .5 
DATA  FPS/ .000005/ 

DATA  PIF/3. 1415927/ 

PHT(X)=1. /SQ PT( 2  *P I E ) *EXP( —X*  *  2/2 . ) 

P(  X  ,C. )  =(  1  .+r)/(  (X* *2  +  2. /PTE*  (  1  .+G)**2)  **  .5  +  G*X) 

PAM I N( X ) - . 7 1  *  X/* 1 . +  .7 1 *X ) 

TPL  (  X  ,GI,  ,f  >U )  =  (  PAMIN  (  X)*P(X,GL)  +  (  1  .-PAMIN(X)  )*P(  X,GU)  )*PHI  (X) 
GMTN-PI E- 1 . 

GMAX=2 . / ( PIE-2 . ) 

X=1  . 

102  TP=TPL( X ,GMIN ,GKAX ) 

IF  (TP.LT.  ALPHA)  GO  TO  103 

X=X+1 

GO  TO  102 

103  Rlr—  X—  1  . 

BU=X 

104  X=  (  BL+BU )  /'2  . 

TP=TPL( X ,GMIN ,GMAX ) 

IF  (TP.LT.  ALPHA)  BU=X 
IF  (TP.GE.  ALPHA)  BL=X 
IF  (ABS(BU-BL) . LE.  EPS)  GO  TO  105 
GO  TO  104 

105  XALP=( BL+BU )/2. 

RETURN 

END 


SUBROUTINE"  TT(DF,X,TP) 

C  T  TAIL  AREA  FROM  X  TO  INFINITY,  X  NON- NEGATIVE 

P(X,C,n)=( 1 .+C)/( (X**2+4.*CK(P) **2* ( 1 .+G)**2)**.5+G*X) 
FKX(X,P)=CK(P)*(  1 . +  X**  2/P ) ** ( (-P+1 . )/2. ) 

GMAX ( X ) =4 . *CK ( X ) *  *2/ ( 1 . -4 . *CK ( X ) **2 ) 

GMIN ( X )=X/ ( 2 . * ( X+2 . ) *CK( X ) **2 )- 1 . 

TCON (  D  F )  = .  7  1  +  1 .8/DF+1 .9/DF**2 
PAMIN(DF,X)=TCON(DF)*X/( 1 . +TCON ( PF ) *X ) 

TPL(DF,X,GL,GU)=(PAMIN(DF,X)*P(X,GL,DF)+( 1 PAMIN(DF,X) )* 

1 P( X ,GU , DF ) ) *FKX ( X , PF ) 

TP=TPL( DF ,X ,GMIN ( DF ) ,GMAX( DF) ) 

RETURN 

FNP 

SUBROUTINE  TTI  (  PF,  AI.PHA  ,X  ALP) 

C  UPPER  T  PERCENTILE  XALP  OORRF.SPONP1NG  TO  TAIL  PROBABILITY  ALPHA, 
C  ALPHA  BETWEEN  0  AND  .5 

P(X,G,D)  =  (  1  .+G)/(  ( X*  *2+4 .  *CK  {  D)  **2*  (  1  . +G)  **2 )  *  * .  5+C-*X ) 

FKX (X,P)=CK(D)*( 1 .+X**2/D) ** ( (-D+1. )/2. ) 

GMAX ( X ) =4 . *CK( X ) **2/ ( 1 .-4 . *OK ( X ) * *2 ) 

GMTN(X)=X/(2.*(X+2. ) *CX { X ) **2 ) - 1 . 

TCON ( PF) = . 7 1+ 1 . R/PF+ 1 . 9/PF**2 

PAMT  N ( PF , X )  TCON(PF)*X/(  1  . +TCON ( PF ) *X ) 

TPL(  PF ,  X  ,GI.  ,GU  )  -  (  PAMI N  (PE",X)*P(X,GL,PF)  +  (  1  PAMIN(  DF,X)  )* 

1 P( X ,GU ,PF ) ) *FKX ( X ,DF ) 

PATA  EPS/. 000005/ 

X=1  . 

102  TP=TPL( DF,X ,GMIN( PF) ,GMAX(PF) ) 

IF  (TP.LT.  ALPHA)  GO  TO  103 
X=X+1 . 

GO  TO  102 

103  BL=X- 1 . 

BU=X 

104  X  -  (  BI.+  BU  )  / 2  . 

TP-TPI,(  PF , X  ,GMTN(  PF)  , C. MAX  (  PF)  ) 

IF  (TP.LT.  ALPHA)  BU=X 
IF  (TP.GE.  ALPHA)  BL=X 
I F ( AES ( BU-BL ) . LE .  EPS)  GO  TO  105 
GO  TO  104 

105  XALP= ( BL+BU ) /2 . 

RETURN 

END 

FUNCTION  CK(X) 

ALNG ( X ) = ( X- . 5 ) *ALOG ( X ) -X+ 1/ (  1 2 . *X ) -  1/ ( 360 . *X* *3 ) 

GAMMA (T)=1 . /T-. 577 19 165+ . 988 20589*T-. 8970 5694*T**2+.9 1820686* 
IT* *3-. 75670408* T**4+. 48219939* T* *5-. 19352782*T**6+. 03586834* 
2T**7 

PAT,,  PTE/3 .141  5927/ 

Zl=(X+1)/2. 

Z2-X/2 . 


»  n  *v  ■!  » ;u  ■  v  ■  u  ■  i  »k 


■i  ■  i»  ■  i 


i ■  "i  •  "i  >»  '  «  p» 


rrm 


l  ■  L 


IF  (Z2  .GT.  10. )  GO  TO  11 
TEMP 1=1. 

201  IF  (Zl  .LE.  1.)  GO  TO  202 
TEMP 1 =TEMP 1*(Z1— 1. ) 

Z1=Z1- 1 . 

GO  TO  201 

202  TEMP1=TEMP1*GAMMA(Z1 ) 

TEMP2= 1 . 

203  IF  (Z2  .LE.  1 . )  GO  TO  204 
TEMP2=TEMP2* ( Z2- 1 . ) 

Z2=Z2- 1 . 

GO  TO  203 

204  TEMP2=TEMP2*GAMMA(Z2) 

CK=TEMP 1 /TEMP 2 /SQ RT ( P I E* X ) 

RETURN 

11  CK=EXP( ALNG( ( X+ 1 . )/2 . ) -ALNG( X/2 . ) - . 5*ALOG( PIE*X ) ) 
RETURN 
END 
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